Now that we have way too many copeopod time series, lets see if any of them are useful covariates in the Atlantic herring assessment model.
The primary target is a covariate on Atlantic herring recruitment.
Adelle’s boosted regression tree model shows both spring large copepods and fall small copepods as important factors explaining recruitment patterns estimated by the WHAM model.
We have also developed a herring-larvae specific small copepod index within the area with more than the 70th percentile of cumulative 1982-2022 herring larvae density.
So we can try at least three flavors of copepod index as a model covariate.
The code that takes the output of the VAST models and makes WHAM inputs is WHAMinputs.R
Also adding a zooplankton volume spring index to try in the place of the spring large copepods.
Installed multiWHAM from the devel branch: https://github.com/timjmiller/wham/tree/devel
devtools::install_github("timjmiller/wham", dependencies=TRUE, ref="devel")
Installation worked on the outdated mac os… so far.
Lets test the installation by seeing if I can get the same outputs as Jon did by running the model without any additions.
Modify Jon’s code from the Atlantic Herring RTWG google drive and now
in this directory, Example WHAM code to share.R which
sources his WHAMfxns.R and hindcast.R also now
in this directory.
First run the base model and compare with the stored to ensure my installation is ok:
#source some functions
source(here::here("WHAMfits/WHAMfxns.R"))
#source the hindcast function for MASE etc.
source(here::here("WHAMfits/hindcast.R"))
# starting with mm192, add suffix for different attempts
# test with no change
config <- "nochange"
# name the model directory
name <- paste0("mm192_", config)
write.dir <- here::here(sprintf("WHAMfits/%s/", name))
if(!dir.exists(write.dir)) {
dir.create(write.dir)
}
# copy the input dat file to model directory
asapRun2_87 <- here::here("WHAMfits/mm192/Run2_87.DAT")
file.copy(from=file.path(asapRun2_87), to=write.dir, overwrite=FALSE)
asapRun2dat<- read_asap3_dat(here::here(write.dir, "Run2_87.DAT"))
# need data files read in to run these--but don't,
# the full data object is already available in the mm192 rds file
#asapRun2dat=CVfxn(asapRun9dat) #empirical CVs
#asapRun2dat=ESSfxn(asapRun9dat,addnumcatch=150,addnumsurv=150)
mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))
input <- mm192mod$input
m1 <- fit_wham(input, do.osa = F)
saveRDS(m1, file.path(write.dir, "mm192nochange.rds"))
check_convergence(m1)
# these should be identical; they are
test <- compare_wham_models(mods = list("orig" = mm192mod, "nochange" = m1), fdir = write.dir)
saveRDS(test, file.path(write.dir, "compare.rds"))
The original model and my rerun using the same inputs (nochange) directly overlap:
knitr::include_graphics(here::here("WHAMfits/mm192_nochange/compare_png/compare_SSB_F_R.png"))
New explore environmental covariate (ecov) options. What format does the input need?
Following https://timjmiller.github.io/wham/articles/ex2_CPI_recruitment.html columns Year, Index, Index_sigma
We are testing indices of bottom-up recruitment impacts–food for larvae, postlarvae, and juveniles
Format ecov indices; we’ll try at least 3 variants:
The code in WHAMinputs.R was used to create csv files of the VAST derives zooplankton indices. Files are saved in the WHAMfits folder.
Start with spring large copepods. Following the suggestions in WHAM vignette 2:
OK thats outmoded. Can just add an ecov object to existing input
using set_ecov:
config <- "largecope"
# name the model directory
name <- paste0("mm192_", config)
write.dir <- here::here(sprintf("WHAMfits/%s/", name))
if(!dir.exists(write.dir)) {
dir.create(write.dir)
}
mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))
input <- mm192mod$input
# asapdat<- read_asap3_dat(here::here("WHAMfits/mm192/Run2_87.DAT"))
#
# asap3 <- input$asap3
env.dat <- read.csv(here::here("WHAMfits/springlargecopeindex.csv"), header=T)
# make it the same length as model ?
#env.dat <- env.dat[env.dat$Time > 1986,]
# test a single ecov setup with no effect first
ecov_on <- list(
label = "lgCopeSpr",
mean = as.matrix(env.dat$her_sp_Estimate/1000000),
logsigma = as.matrix(log(env.dat$her_sp_SE/1000000)), #"est_1",
year = env.dat$Time,
use_obs = matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
lag = 0, # postlarval and juvenile food in year t affects recruitment in year t
process_model = "ar1", # "rw" or "ar1"
#where = "recruit",
recruitment_how = as.matrix("controlling-lag-0-linear")) # 0 = no effect, 1 = controlling, 2 = limiting
# modify current model input call adding ecov list
#input$call
ecovinput <- set_ecov(input, ecov=ecov_on)
# ecovinput <- prepare_wham_input(asap3 = asapdat, # correct structurr, asap3 is wrong structure
# model_name = "test", #df.mods[m, "Model"],
# ecov = ecov,
# recruit_model = 2,
# selectivity = list(fix_pars = list(7:8,
# 2,
# c(1:2, 4:8),
# 1:8,
# c(1, 3:8),
# NULL,
# c(1, 4:8),
# NULL)),
# M = list(re_model = as.matrix("none")), #as.matrix(df.mods[m, "M_re"])),
# NAA_re = list(sigma = "rec+1",
# cor = "ar1_a", #df.mods[m, "naa_cor"],
# N1_model = "age-specific-fe", #init_Nnum[m],
# decouple_recruitment = TRUE),
# age_comp = "logistic-normal-ar1-miss0",
# basic_info = list(bias_correct_process = FALSE,
# bias_correct_BRPs = FALSE,
# percentSPR = 40,
# percentFXSPR = 100,
# XSPR_R_avg_yrs = 6:35,
# XSPR_R_opt = 5))
ecovon_mod <- fit_wham(ecovinput, do.osa = F)
saveRDS(ecovon_mod, file.path(write.dir, "testecovon.rds"))
plot_wham_output(mod = ecovon_mod, dir.main = file.path(write.dir))
# test a single ecov setup with no effect first
ecov_off <- list(
label = "lgCopeSpr",
mean = as.matrix(env.dat$her_sp_Estimate/1000000),
logsigma = as.matrix(log(env.dat$her_sp_SE/1000000)),
year = env.dat$Time,
use_obs = matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
lag = 0, # postlarval and juvenile food in year t affects recruitment in year t
process_model = "ar1", # "rw" or "ar1"
#where = "recruit",
recruitment_how = as.matrix("none")) # 0 = no effect, 1 = controlling, 2 = limiting
# modify current model input call adding ecov list
#input$call
ecovinput <- set_ecov(input, ecov=ecov_off)
mod <- fit_wham(ecovinput, do.osa = F)
saveRDS(mod, file.path(write.dir, "testecovoff.rds"))
#compare to original, not the same data
test <- compare_wham_models(mods = list("testecovoff" = mod, "testecovon" = ecovon_mod), fdir = write.dir)
saveRDS(test, file.path(write.dir, "testcompare.rds"))
Test run was broken until I changed units of zooplankton to millions of cells.
Now everything runs, the covariate is included, and makes no difference in the model (which looks identical to the base model so thats a plus):
testcompare <- readRDS(here::here("WHAMfits/mm192_largecope/testcompare.rds"))
knitr::include_graphics(here::here("WHAMfits/mm192_largecope/plots_png/results/Ecov_1_lgCopeSpr.png"))
testcompare$tab
## dAIC AIC rho_R rho_SSB rho_Fbar
## testecovoff 0.0 -2759.1 0.8682 0.5901 -0.2314
## testecovon 0.4 -2758.7 0.9081 0.5920 -0.2379
testcompare$g[[1]]
Lets try a suite of models with the same recruitment process but different ways the large copepod series relates to it.
We’ll always use recruitment process 2, random about mean, because that is in mm192.
We’ll look at both random walk “rw” and AR1 “ar1” processes for the environmental covariate.
We’ll set the parameterrecruitment_how as one of the
following:
{= “none”}{no effect.}
{= “controlling”}{pre-recruit density-independent mortality.}
{= “limiting”}{ maximum recruitment, e.g. ecov determines amount of
suitable habitat)}
{= “lethal”}{threshold, i.e. R –> 0 at some ecov value.}
{= “masking”}{metabolic/growth, decreases dR/dS}
{= “directive”}{e.g. behavioral}
Going to try “none”, “controlling”, “limiting”, and “masking” (note that “masking” didn’t work well in the state space RTA simulations).
config <- "largecope"
# name the model directory
name <- paste0("mm192_", config)
write.dir <- here::here(sprintf("WHAMfits/%s/", name))
if(!dir.exists(write.dir)) {
dir.create(write.dir)
}
mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))
input <- mm192mod$input
# asapdat<- read_asap3_dat(here::here("WHAMfits/mm192/Run2_87.DAT"))
#
# asap3 <- input$asap3
env.dat <- read.csv(here::here("WHAMfits/springlargecopeindex.csv"), header=T)
# larger set of ecov setups to compare
df.mods <- data.frame(Recruitment = c(rep(c(2,2,3,3), 2)),
ecov_process = c(rep("rw",4),rep("ar1",4)),
ecov_how = rep(c("none",
"controlling-lag-0-linear",
"limiting-lag-0-linear",
"masking-lag-0-linear"),2), stringsAsFactors=FALSE)
n.mods <- dim(df.mods)[1]
df.mods$Model <- paste0("m",1:n.mods)
df.mods <- dplyr::select(df.mods, Model, tidyselect::everything()) # moves Model to first col
for(m in 5:n.mods){
ecov <- list(
label = "lgCopeSpr",
mean = as.matrix(env.dat$her_sp_Estimate/1000000),
logsigma = as.matrix(log(env.dat$her_sp_SE/1000000)),
year = env.dat$Time,
use_obs = matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
#lag = 0, # postlarval and juvenile food in year t affects recruitment in year t
process_model = df.mods$ecov_process[m], # "rw" or "ar1"
#where = c("recruit")[as.logical(df.mods$ecov_how[m])+1],
recruitment_how = as.matrix(df.mods$ecov_how[m])
) # 0 = no effect, 1 = controlling, 2 = limiting
input$data$recruit_model <- df.mods$Recruitment[m]
input$options$basic_info$recruit_model <- df.mods$Recruitment[m]
ecovinput <- set_ecov(input, ecov=ecov)
mod <- fit_wham(ecovinput, do.osa = T)
# Save model
saveRDS(mod, file.path(write.dir, paste0(df.mods$Model[m],".rds")))
# Plot output in new subfolder
plot_wham_output(mod=mod, dir.main=file.path(write.dir,df.mods$Model[m]), out.type='html')
}
Check all for convergence (from https://timjmiller.github.io/wham/articles/ex2_CPI_recruitment.html)
config <- "largecope"
# name the model directory
name <- paste0("mm192_", config)
write.dir <- here::here(sprintf("WHAMfits/%s/", name))
# need this for html to work because block above running model not evaluated on knit
df.mods <- data.frame(Recruitment = c(rep(c(2,2,3,3), 2)),
ecov_process = c(rep("rw",4),rep("ar1",4)),
recruitment_how = rep(c("none",
"controlling-lag-0-linear",
"limiting-lag-0-linear",
"masking-lag-0-linear"),2), stringsAsFactors=FALSE)
n.mods <- dim(df.mods)[1]
df.mods$Model <- paste0("m",1:n.mods)
df.mods <- dplyr::select(df.mods, Model, tidyselect::everything()) # moves Model to first col
mod.list <- paste0(write.dir, df.mods$Model,".rds")
mods <- lapply(mod.list, readRDS)
n.mods <- length(mod.list)
vign2_conv <- lapply(mods, function(x) capture.output(check_convergence(x)))
for(m in 1:n.mods) cat(paste0("Model ",m,":"), vign2_conv[[m]], "", sep='\n')
## Model 1:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 7.31e-12
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 2:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 4.47e-03
## Max gradient parameter: logit_q
## TMB:sdreport() was performed successfully for this model
##
## Model 3:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 4.96e-03
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 4:
## stats:nlminb thinks the model has NOT converged: mod$opt$convergence != 0
## Maximum gradient component: 5.38e+07
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 5:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 9.33e-12
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 6:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.37e-03
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed successfully for this model
##
## Model 7:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.88e-02
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 8:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.62e-03
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed successfully for this model
Compare models (from https://timjmiller.github.io/wham/articles/ex2_CPI_recruitment.html)
opt_conv = 1-sapply(mods, function(x) x$opt$convergence)
ok_sdrep = sapply(mods, function(x) if(x$na_sdrep==FALSE & !is.na(x$na_sdrep)) 1 else 0)
df.mods$conv <- as.logical(opt_conv)
df.mods$pdHess <- as.logical(ok_sdrep)
df.mods$NLL <- sapply(mods, function(x) round(x$opt$objective,3))
not_conv <- !df.mods$conv | !df.mods$pdHess
mods2 <- mods
mods2[not_conv] <- NULL
df.aic.tmp <- as.data.frame(compare_wham_models(mods2, table.opts=list(sort=FALSE, calc.rho=T))$tab)
## dAIC AIC rho_R rho_SSB rho_Fbar
## m1 3.7 -2755.4 0.8594 0.5907 -0.2362
## m2 4.3 -2754.8 0.8104 0.5940 -0.2448
## m3 4.3 -2754.8 0.8162 0.6011 -0.2446
## m4 0.0 -2759.1 0.8682 0.5901 -0.2314
## m5 0.4 -2758.7 0.9081 0.5920 -0.2379
## m6 0.4 -2758.7 0.5715 0.5914 -0.2399
## m7 1.9 -2757.2 0.8908 0.5945 -0.2374
df.aic <- df.aic.tmp[FALSE,]
ct = 1
for(i in 1:n.mods){
if(not_conv[i]){
df.aic[i,] <- rep(NA,5)
} else {
df.aic[i,] <- df.aic.tmp[ct,]
ct <- ct + 1
}
}
df.mods <- cbind(df.mods, df.aic)
df.mods <- df.mods[order(df.mods$dAIC, na.last=TRUE),]
df.mods[is.na(df.mods$AIC), c('dAIC','AIC','rho_R','rho_SSB','rho_Fbar')] <- "---"
rownames(df.mods) <- NULL
df.mods
## Model Recruitment ecov_process recruitment_how conv pdHess
## 1 m5 2 ar1 none TRUE TRUE
## 2 m6 2 ar1 controlling-lag-0-linear TRUE TRUE
## 3 m7 3 ar1 limiting-lag-0-linear TRUE TRUE
## 4 m8 3 ar1 masking-lag-0-linear TRUE TRUE
## 5 m1 2 rw none TRUE TRUE
## 6 m2 2 rw controlling-lag-0-linear TRUE TRUE
## 7 m3 3 rw limiting-lag-0-linear TRUE TRUE
## 8 m4 3 rw masking-lag-0-linear FALSE FALSE
## NLL dAIC AIC rho_R rho_SSB rho_Fbar
## 1 -1509.570 0 -2759.1 0.8682 0.5901 -0.2314
## 2 -1510.349 0.4 -2758.7 0.9081 0.592 -0.2379
## 3 -1510.349 0.4 -2758.7 0.5715 0.5914 -0.2399
## 4 -1509.592 1.9 -2757.2 0.8908 0.5945 -0.2374
## 5 -1506.709 3.7 -2755.4 0.8594 0.5907 -0.2362
## 6 -1507.424 4.3 -2754.8 0.8104 0.594 -0.2448
## 7 -1507.424 4.3 -2754.8 0.8162 0.6011 -0.2446
## 8 -1175.320 --- --- --- --- ---
Outputs are nearly identical for SSB and F, but there are small differences in status determination.
Ignore results from the limiting and masking models as Bev-Holt parameters aren’t really being estimated.
The numbers in the plots are not the model numbers! They are in order by the first table above
lgcope1 <- compare_wham_models(mods2, do.table= FALSE, fdir = write.dir)
saveRDS(lgcope1, file.path(write.dir, "lgcope1compare.rds"))
lgcope1$g[[1]]
lgcope1$g[[8]]
lgcope1$g[[9]]
knitr::include_graphics(here::here(write.dir, "compare_png/compare_rel_status_kobe.png"))
Similar approach for first pass at fall small copepods. These will lag 1.
Setup
config <- "smcopeFall"
# name the model directory
name <- paste0("mm192_", config)
write.dir <- here::here(sprintf("WHAMfits/%s/", name))
if(!dir.exists(write.dir)) {
dir.create(write.dir)
}
mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))
input <- mm192mod$input
# asapdat<- read_asap3_dat(here::here("WHAMfits/mm192/Run2_87.DAT"))
#
# asap3 <- input$asap3
# larger set of ecov setups to compare
df.mods <- data.frame(Recruitment = c(rep(c(2,2,3,3), 2)),
ecov_process = c(rep("rw",4),rep("ar1",4)),
ecov_how = rep(c("none",
"controlling-lag-1-linear",
"limiting-lag-1-linear",
"masking-lag-1-linear"),2), stringsAsFactors=FALSE)
n.mods <- dim(df.mods)[1]
df.mods$Model <- paste0("m",1:n.mods)
df.mods <- dplyr::select(df.mods, Model, tidyselect::everything()) # moves Model to first col
Run models
env.dat <- read.csv(here::here("WHAMfits/fallsmallcopeALLindex.csv"), header=T)
# 2020 is missing
env.dat[env.dat == 0] <- NA
# don't use the NA value
use.obs <- matrix(1, ncol=1, nrow=dim(env.dat)[1])
use.obs[39,] <- 0
for(m in 1:n.mods){
ecov <- list(
label = "msCopeFall",
mean = as.matrix(env.dat$her_fa_Estimate/1000000),
logsigma = as.matrix(log(env.dat$her_fa_SE/1000000)),
year = env.dat$Time,
use_obs = use.obs, # matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
#lag = 0, # postlarval and juvenile food in year t affects recruitment in year t
process_model = df.mods$ecov_process[m], # "rw" or "ar1"
#where = c("recruit")[as.logical(df.mods$ecov_how[m])+1],
recruitment_how = as.matrix(df.mods$ecov_how[m])
) # 0 = no effect, 1 = controlling, 2 = limiting
# this isn't working, not estimating a recruitment function
# for the best because it will break
# for reference, this code should be added to attempt
#input206 <- set_NAA(mm192$mm192$input, NAA_re=list(recruit_model=3,sigma="rec+1",cor="ar1_a",N1_model="age-specific-fe",decouple_recruitment=TRUE))
input$data$recruit_model <- df.mods$Recruitment[m]
input$options$basic_info$recruit_model <- df.mods$Recruitment[m]
ecovinput <- set_ecov(input, ecov=ecov)
mod <- fit_wham(ecovinput, do.osa = T)
# Save model
saveRDS(mod, file.path(write.dir, paste0(df.mods$Model[m],".rds")))
# Plot output in new subfolder
plot_wham_output(mod=mod, dir.main=file.path(write.dir,df.mods$Model[m]), out.type='html')
}
Estimated small copepods all Fall covariate (“ar1”, m5):
knitr::include_graphics(here::here(write.dir, "m5/plots_png/results/Ecov_1_msCopeFall.png"))
Compare models (from https://timjmiller.github.io/wham/articles/ex2_CPI_recruitment.html)
Ignore results from the limiting and masking models as Bev-Holt parameters aren’t really being estimated.
Also nothing worked that had a recruitment effect
mod.list <- paste0(write.dir, df.mods$Model,".rds")
mods <- lapply(mod.list, readRDS)
n.mods <- length(mod.list)
vign2_conv <- lapply(mods, function(x) capture.output(check_convergence(x)))
for(m in 1:n.mods) cat(paste0("Model ",m,":"), vign2_conv[[m]], "", sep='\n')
## Model 1:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.20e-12
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 2:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 4.65e-01
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 3:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.06e+00
## Max gradient parameter: F_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 4:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.34e-08
## Max gradient parameter: F_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 5:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.39e-11
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 6:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 7.34e+00
## Max gradient parameter: logit_selpars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 7:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.64e-03
## Max gradient parameter: F_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 8:
## stats:nlminb thinks the model has NOT converged: mod$opt$convergence != 0
## Maximum gradient component: 1.15e+07
## Max gradient parameter: F_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
opt_conv = 1-sapply(mods, function(x) x$opt$convergence)
ok_sdrep = sapply(mods, function(x) if(x$na_sdrep==FALSE & !is.na(x$na_sdrep)) 1 else 0)
df.mods$conv <- as.logical(opt_conv)
df.mods$pdHess <- as.logical(ok_sdrep)
df.mods$NLL <- sapply(mods, function(x) round(x$opt$objective,3))
not_conv <- !df.mods$conv | !df.mods$pdHess
mods2 <- mods
mods2[not_conv] <- NULL
df.aic.tmp <- as.data.frame(compare_wham_models(mods2, table.opts=list(sort=FALSE, calc.rho=T))$tab)
## dAIC AIC rho_R rho_SSB rho_Fbar
## m1 11.6 -2597.8 3.2038 0.6009 -0.2487
## m2 0.0 -2609.4 0.8624 0.5830 -0.2315
df.aic <- df.aic.tmp[FALSE,]
ct = 1
for(i in 1:n.mods){
if(not_conv[i]){
df.aic[i,] <- rep(NA,5)
} else {
df.aic[i,] <- df.aic.tmp[ct,]
ct <- ct + 1
}
}
df.mods <- cbind(df.mods, df.aic)
df.mods <- df.mods[order(df.mods$dAIC, na.last=TRUE),]
df.mods[is.na(df.mods$AIC), c('dAIC','AIC','rho_R','rho_SSB','rho_Fbar')] <- "---"
rownames(df.mods) <- NULL
df.mods
## Model Recruitment ecov_process ecov_how conv pdHess
## 1 m5 2 ar1 none TRUE TRUE
## 2 m1 2 rw none TRUE TRUE
## 3 m2 2 rw controlling-lag-1-linear TRUE FALSE
## 4 m3 3 rw limiting-lag-1-linear TRUE FALSE
## 5 m4 3 rw masking-lag-1-linear TRUE FALSE
## 6 m6 2 ar1 controlling-lag-1-linear TRUE FALSE
## 7 m7 3 ar1 limiting-lag-1-linear TRUE FALSE
## 8 m8 3 ar1 masking-lag-1-linear FALSE FALSE
## NLL dAIC AIC rho_R rho_SSB rho_Fbar
## 1 -1434.708 0 -2609.4 0.8624 0.583 -0.2315
## 2 -1427.902 11.6 -2597.8 3.2038 0.6009 -0.2487
## 3 -1429.056 --- --- --- --- ---
## 4 -1429.056 --- --- --- --- ---
## 5 -1427.902 --- --- --- --- ---
## 6 -1426.200 --- --- --- --- ---
## 7 -1428.500 --- --- --- --- ---
## 8 -1431.880 --- --- --- --- ---
smcope1 <- compare_wham_models(mods2, do.table= FALSE, fdir = write.dir)
saveRDS(smcope1, file.path(write.dir, "smcope1compare.rds"))
smcope1$g[[1]]
smcope1$g[[8]]
smcope1$g[[9]]
knitr::include_graphics(here::here(write.dir, "compare_png/compare_rel_status_kobe.png"))
Similar approach for first pass at fall small copepods. These will lag 1.
Setup
config <- "smcopeSepFeb"
# name the model directory
name <- paste0("mm192_", config)
write.dir <- here::here(sprintf("WHAMfits/%s/", name))
if(!dir.exists(write.dir)) {
dir.create(write.dir)
}
mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))
input <- mm192mod$input
# asapdat<- read_asap3_dat(here::here("WHAMfits/mm192/Run2_87.DAT"))
#
# asap3 <- input$asap3
# larger set of ecov setups to compare
df.mods <- data.frame(Recruitment = c(rep(c(2,2,3,3), 2)),
ecov_process = c(rep("rw",4),rep("ar1",4)),
ecov_how = rep(c("none",
"controlling-lag-1-linear",
"limiting-lag-1-linear",
"masking-lag-1-linear"),2), stringsAsFactors=FALSE)
n.mods <- dim(df.mods)[1]
df.mods$Model <- paste0("m",1:n.mods)
df.mods <- dplyr::select(df.mods, Model, tidyselect::everything()) # moves Model to first col
Run models
env.dat <- read.csv(here::here("WHAMfits/sepfebsmallcopeALLlarvareaindex.csv"), header=T)
# 2020 is missing
env.dat[env.dat == 0] <- NA
# don't use the NA value
use.obs <- matrix(1, ncol=1, nrow=dim(env.dat)[1])
use.obs[39,] <- 0
for(m in 1:n.mods){
ecov <- list(
label = "smCopeSepFeb",
mean = as.matrix(env.dat$her_larv_Estimate/1000000),
logsigma = as.matrix(log(env.dat$her_larv_SE/1000000)),
year = env.dat$Time,
use_obs = use.obs, # matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
#lag = 0, # postlarval and juvenile food in year t affects recruitment in year t
process_model = df.mods$ecov_process[m], # "rw" or "ar1"
#where = c("recruit")[as.logical(df.mods$ecov_how[m])+1],
recruitment_how = as.matrix(df.mods$ecov_how[m])
) # 0 = no effect, 1 = controlling, 2 = limiting
# this isn't working, not estimating a recruitment function
# for the best because it will break
# for reference, this code should be added to attempt
#input206 <- set_NAA(mm192$mm192$input, NAA_re=list(recruit_model=3,sigma="rec+1",cor="ar1_a",N1_model="age-specific-fe",decouple_recruitment=TRUE))
input$data$recruit_model <- df.mods$Recruitment[m]
input$options$basic_info$recruit_model <- df.mods$Recruitment[m]
ecovinput <- set_ecov(input, ecov=ecov)
mod <- fit_wham(ecovinput, do.osa = T)
# Save model
saveRDS(mod, file.path(write.dir, paste0(df.mods$Model[m],".rds")))
# Plot output in new subfolder
plot_wham_output(mod=mod, dir.main=file.path(write.dir,df.mods$Model[m]), out.type='html')
}
Estimated small copepods September-February in herring larval area covariate (“ar1”, m5):
knitr::include_graphics(here::here(write.dir, "m5/plots_png/results/Ecov_1_smCopeSepFeb.png"))
Compare models (from https://timjmiller.github.io/wham/articles/ex2_CPI_recruitment.html)
Ignore results from the limiting and masking models as Bev-Holt parameters aren’t really being estimated.
Also nothing worked that had a recruitment effect
mod.list <- paste0(write.dir, df.mods$Model,".rds")
mods <- lapply(mod.list, readRDS)
n.mods <- length(mod.list)
vign2_conv <- lapply(mods, function(x) capture.output(check_convergence(x)))
for(m in 1:n.mods) cat(paste0("Model ",m,":"), vign2_conv[[m]], "", sep='\n')
## Model 1:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.88e-11
## Max gradient parameter: F_pars
## TMB:sdreport() was performed successfully for this model
##
## Model 2:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.51e+00
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 3:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.21e+00
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 4:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 5.99e+00
## Max gradient parameter: Ecov_process_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 5:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 7.46e-12
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 6:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.37e+00
## Max gradient parameter: F_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 7:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 5.91e-04
## Max gradient parameter: logit_q
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 8:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 9.60e-09
## Max gradient parameter: F_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
opt_conv = 1-sapply(mods, function(x) x$opt$convergence)
ok_sdrep = sapply(mods, function(x) if(x$na_sdrep==FALSE & !is.na(x$na_sdrep)) 1 else 0)
df.mods$conv <- as.logical(opt_conv)
df.mods$pdHess <- as.logical(ok_sdrep)
df.mods$NLL <- sapply(mods, function(x) round(x$opt$objective,3))
not_conv <- !df.mods$conv | !df.mods$pdHess
mods2 <- mods
mods2[not_conv] <- NULL
df.aic.tmp <- as.data.frame(compare_wham_models(mods2, table.opts=list(sort=FALSE, calc.rho=T))$tab)
## dAIC AIC rho_R rho_SSB rho_Fbar
## m1 3.9 -2627.3 0.8541 0.5879 -0.2311
## m2 0.0 -2631.2 0.8793 0.6035 -0.2472
df.aic <- df.aic.tmp[FALSE,]
ct = 1
for(i in 1:n.mods){
if(not_conv[i]){
df.aic[i,] <- rep(NA,5)
} else {
df.aic[i,] <- df.aic.tmp[ct,]
ct <- ct + 1
}
}
df.mods <- cbind(df.mods, df.aic)
df.mods <- df.mods[order(df.mods$dAIC, na.last=TRUE),]
df.mods[is.na(df.mods$AIC), c('dAIC','AIC','rho_R','rho_SSB','rho_Fbar')] <- "---"
rownames(df.mods) <- NULL
df.mods
## Model Recruitment ecov_process ecov_how conv pdHess NLL
## 1 m5 2 ar1 none TRUE TRUE -1445.582
## 2 m1 2 rw none TRUE TRUE -1442.652
## 3 m2 2 rw controlling-lag-1-linear TRUE FALSE -1444.922
## 4 m3 3 rw limiting-lag-1-linear TRUE FALSE -1444.922
## 5 m4 3 rw masking-lag-1-linear TRUE FALSE -1444.215
## 6 m6 2 ar1 controlling-lag-1-linear TRUE FALSE -1447.832
## 7 m7 3 ar1 limiting-lag-1-linear TRUE FALSE -1439.374
## 8 m8 3 ar1 masking-lag-1-linear TRUE FALSE -1445.582
## dAIC AIC rho_R rho_SSB rho_Fbar
## 1 0 -2631.2 0.8793 0.6035 -0.2472
## 2 3.9 -2627.3 0.8541 0.5879 -0.2311
## 3 --- --- --- --- ---
## 4 --- --- --- --- ---
## 5 --- --- --- --- ---
## 6 --- --- --- --- ---
## 7 --- --- --- --- ---
## 8 --- --- --- --- ---
smcopelarv1 <- compare_wham_models(mods2, do.table= FALSE, fdir = write.dir)
saveRDS(smcopelarv1, file.path(write.dir, "smcopelarv1compare.rds"))
smcopelarv1$g[[1]]
smcopelarv1$g[[8]]
smcopelarv1$g[[9]]
knitr::include_graphics(here::here(write.dir, "compare_png/compare_rel_status_kobe.png"))
Try the spring zooplankton with full zooplankton volume instead of just large copepods.
Setup
config <- "zoopvolSpring"
# name the model directory
name <- paste0("mm192_", config)
write.dir <- here::here(sprintf("WHAMfits/%s/", name))
if(!dir.exists(write.dir)) {
dir.create(write.dir)
}
mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))
input <- mm192mod$input
# asapdat<- read_asap3_dat(here::here("WHAMfits/mm192/Run2_87.DAT"))
#
# asap3 <- input$asap3
# larger set of ecov setups to compare
df.mods <- data.frame(Recruitment = c(rep(c(2,2,3,3), 2)),
ecov_process = c(rep("rw",4),rep("ar1",4)),
ecov_how = rep(c("none",
"controlling-lag-0-linear",
"limiting-lag-0-linear",
"masking-lag-0-linear"),2), stringsAsFactors=FALSE)
n.mods <- dim(df.mods)[1]
df.mods$Model <- paste0("m",1:n.mods)
df.mods <- dplyr::select(df.mods, Model, tidyselect::everything()) # moves Model to first col
Run models
env.dat <- read.csv(here::here("WHAMfits/springzoopvolumeindex.csv"), header=T)
# units of zoop volume are ml so lets convert to liters?
for(m in 1:n.mods){
ecov <- list(
label = "zoopvolSpr",
mean = as.matrix(env.dat$her_sp_Estimate/1000),
logsigma = as.matrix(log(env.dat$her_sp_SE/1000)),
year = env.dat$Time,
use_obs = matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
#lag = 0, # postlarval and juvenile food in year t affects recruitment in year t
process_model = df.mods$ecov_process[m], # "rw" or "ar1"
#where = c("recruit")[as.logical(df.mods$ecov_how[m])+1],
recruitment_how = as.matrix(df.mods$ecov_how[m])
) # 0 = no effect, 1 = controlling, 2 = limiting
# this isn't working, not estimating a recruitment function
# for the best because it will break
# for reference, this code should be added to attempt
#input <- set_NAA(mm192$mm192$input, NAA_re=list(recruit_model=3,sigma="rec+1",cor="ar1_a",N1_model="age-specific-fe",decouple_recruitment=TRUE))
input$data$recruit_model <- df.mods$Recruitment[m]
input$options$basic_info$recruit_model <- df.mods$Recruitment[m]
ecovinput <- set_ecov(input, ecov=ecov)
mod <- fit_wham(ecovinput, do.osa = T)
# Save model
saveRDS(mod, file.path(write.dir, paste0(df.mods$Model[m],".rds")))
# Plot output in new subfolder
plot_wham_output(mod=mod, dir.main=file.path(write.dir,df.mods$Model[m]), out.type='html')
}
Estimated zooplankton volume covariate (“ar1”, m5):
knitr::include_graphics(here::here(write.dir, "m5/plots_png/results/Ecov_1_zoopvolSpr.png"))
Compare models (from https://timjmiller.github.io/wham/articles/ex2_CPI_recruitment.html)
Ignore results from the limiting and masking models as Bev-Holt parameters aren’t really being estimated.
Also nothing worked that had a recruitment effect
# only take those that are there, but doesnt work with rest of code
#mod.list <- paste0(list.dirs(write.dir, recursive = FALSE), ".rds")
mod.list <- paste0(write.dir, df.mods$Model,".rds")
mods <- lapply(mod.list, readRDS)
n.mods <- length(mod.list)
vign2_conv <- lapply(mods, function(x) capture.output(check_convergence(x)))
for(m in 1:n.mods) cat(paste0("Model ",m,":"), vign2_conv[[m]], "", sep='\n')
## Model 1:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.40e-11
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 2:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 3.20e+00
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 3:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.41e+00
## Max gradient parameter: F_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 4:
## stats:nlminb thinks the model has NOT converged: mod$opt$convergence != 0
## Maximum gradient component: 5.52e+02
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 5:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.12e-11
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 6:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.14e+00
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 7:
## stats:nlminb thinks the model has NOT converged: mod$opt$convergence != 0
## Maximum gradient component: 5.52e+02
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 8:
## stats:nlminb thinks the model has NOT converged: mod$opt$convergence != 0
## Maximum gradient component: 5.52e+02
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
opt_conv = 1-sapply(mods, function(x) x$opt$convergence)
ok_sdrep = sapply(mods, function(x) if(x$na_sdrep==FALSE & !is.na(x$na_sdrep)) 1 else 0)
df.mods$conv <- as.logical(opt_conv)
df.mods$pdHess <- as.logical(ok_sdrep)
df.mods$NLL <- sapply(mods, function(x) round(x$opt$objective,3))
not_conv <- !df.mods$conv | !df.mods$pdHess
mods2 <- mods
mods2[not_conv] <- NULL
df.aic.tmp <- as.data.frame(compare_wham_models(mods2, table.opts=list(sort=FALSE, calc.rho=T))$tab)
## dAIC AIC rho_R rho_SSB rho_Fbar
## m1 18.2 -2590.5 0.7097 0.7253 -0.3005
## m2 0.0 -2608.7 0.8652 0.5978 -0.2361
df.aic <- df.aic.tmp[FALSE,]
ct = 1
for(i in 1:n.mods){
if(not_conv[i]){
df.aic[i,] <- rep(NA,5)
} else {
df.aic[i,] <- df.aic.tmp[ct,]
ct <- ct + 1
}
}
df.mods <- cbind(df.mods, df.aic)
df.mods <- df.mods[order(df.mods$dAIC, na.last=TRUE),]
df.mods[is.na(df.mods$AIC), c('dAIC','AIC','rho_R','rho_SSB','rho_Fbar')] <- "---"
rownames(df.mods) <- NULL
df.mods
## Model Recruitment ecov_process ecov_how conv pdHess
## 1 m5 2 ar1 none TRUE TRUE
## 2 m1 2 rw none TRUE TRUE
## 3 m2 2 rw controlling-lag-0-linear TRUE FALSE
## 4 m3 3 rw limiting-lag-0-linear TRUE FALSE
## 5 m4 3 rw masking-lag-0-linear FALSE FALSE
## 6 m6 2 ar1 controlling-lag-0-linear TRUE FALSE
## 7 m7 3 ar1 limiting-lag-0-linear FALSE FALSE
## 8 m8 3 ar1 masking-lag-0-linear FALSE FALSE
## NLL dAIC AIC rho_R rho_SSB rho_Fbar
## 1 -1434.361 0 -2608.7 0.8652 0.5978 -0.2361
## 2 -1424.248 18.2 -2590.5 0.7097 0.7253 -0.3005
## 3 -1424.585 --- --- --- --- ---
## 4 -1424.585 --- --- --- --- ---
## 5 -1146.241 --- --- --- --- ---
## 6 -1434.730 --- --- --- --- ---
## 7 -1146.241 --- --- --- --- ---
## 8 -1146.241 --- --- --- --- ---
zoopvol <- compare_wham_models(mods2, do.table= FALSE, fdir = write.dir)
saveRDS(zoopvol, file.path(write.dir, "zoopvolcompare.rds"))
zoopvol$g[[1]]
zoopvol$g[[8]]
zoopvol$g[[9]]
knitr::include_graphics(here::here(write.dir, "compare_png/compare_rel_status_kobe.png"))